#####
# CBP: Display map of border region and map of pre-post change in immigration
#####

rm(list = ls())
library(here)
library(qs)
library(data.table)
library(sf)
library(raster)
library(tidyverse)
library(viridis)
library(patchwork)
library(extrafont)


# Preparations ------------------------------------------------------------

# get the crossing border data
cb <- qread(here("data", "cb.qs")); setDT(cb)

# read cantonal borders
canton_geo <- read_sf(here("data", "mappingdata", "timoGrossenbacherData", "g2k15.shp"))

# read country borders 
country_geo <- read_sf(here("data", "mappingdata", "timoGrossenbacherData", "g2l15.shp"))

# read lakes
lake_geo <- read_sf(here("data", "mappingdata", "timoGrossenbacherData", "g2s15.shp"))

# read in raster of relief (this warning can be safely ignored)
relief <- raster(here("data", "mappingdata", "timoGrossenbacherData", "02-relief-ascii.asc")) %>%
  mask(country_geo) %>%
  as("SpatialPixelsDataFrame") %>%
  as.data.frame() %>%
  rename(value = "X02.relief.ascii")

# read productive area. Use the 2019 data but because it has a different projection
# I need to do a bit of work to make it right
muni_geo <- read_sf(here("data", "mappingdata", "K4_polg20190101_vf", "K4polg20190101vf_ch2007Poly.shp"))
muni_geo15 <- read_sf(here("data", "mappingdata", "timoGrossenbacherData", "gde-1-1-15.shp"))
muni_geo <- st_transform(muni_geo, st_crs(muni_geo15))

# clean up
rm(muni_geo15)
gc()

# variable creation -------------------------------------------------------

# create the pre-post dfA measure
cb[ , dfAChange := mean(dfA98[year %in% 2007:2016]) - mean(dfA98[year %in% 1996:1999]), by = bfs19]

# merge the cb with the mapping_data. First make it into a data.table for easy coding
setDT(muni_geo)[cb[year == 2016, .(bfs19, dfAChange, travelMUNmin, BR)], on = c("GDENR" = "bfs19"), `:=`(dfAChange = dfAChange, travelMUNmin = travelMUNmin, BR = BR)]

# create the 4-way variable of border region & travel trime
muni_geo[ , treatmentGroups := factor(
  fcase(
    BR == 0, NA_character_,
    BR == 1 & travelMUNmin < 15, "< 15 minutes",
    BR == 1 & travelMUNmin %between% c(15, 30), "15-30 minutes",
    BR == 1 & travelMUNmin > 30, NA_character_ #AA changed this and the following lines on 2022-05-10 (see old/mapping_travelTim.R for original version)
  ), levels = c("< 15 minutes", "15-30 minutes"))
]

# create a nicely labelled version of the dfA change 
muni_geo[BR ==1 & travelMUNmin <= 30 , #AA added "& travelMUNmin <= 30 " on 2022-07-19
         dfAChangeQuantile := cut(
           dfAChange, breaks = c(-Inf, 1, 10, Inf), include.lowest = T, dig.lab = 3)
]

# turn back into SF object for plotting
muni_geo <- muni_geo[ , .(st_union(geometry), treatmentGroups, dfAChangeQuantile), by = GDENR] %>% sf::st_as_sf()

# plotting ----------------------------------------------------------------

# Prepare the brewer scale
myColorsG <- RColorBrewer::brewer.pal(7, "Greens")[c(7,5,3)]

map <- ggplot(data = muni_geo) +
  geom_raster(data = relief, inherit.aes = F, aes(x = x, y = y, alpha = value)) +
  scale_alpha(name = "", range = c(.6, .0), guide = F) +
  geom_sf(data = country_geo, color = "black", size = 0.05, fill = "transparent") +
  geom_sf(mapping = aes(fill = treatmentGroups, color = factor(is.na(treatmentGroups))), size = 0.025) +
  # scale_size_manual(values = c(.1, 0), guide = F) +
  scale_color_manual(values = c("white", "transparent"), guide = F) +
  scale_fill_manual(
    values = myColorsG,
    # scale_fill_brewer(
    #   palette = "Greens",
    # na.value = "gray80",
    na.translate = F, # doesn't plot the NA or put it in the legend
    name = "Travel time to border crossings",
    labels = c("< 15 min.", "15-30 min.", "30+ min.", "Outside border region"),
    # direction = -1,
    guide = guide_legend(
      override.aes = list(color = "white"),
      direction = "horizontal",
      keyheight = unit(5, units = "pt"),
      keywidth = unit(50, units = "pt"),
      title.position = 'top',
      title.hjust = 0.5,
      label.hjust = 0.5,
      nrow = 1,
      byrow = F,
      label.position = "bottom")
  ) +
  # geom_sf(data = canton_geo, fill = "transparent", color = "gray50", size = 0.5) +
  geom_sf(data = lake_geo, fill = "#D6F1FF", color = "transparent") +
  labs(x = NULL, y = NULL) +
  scale_x_continuous(expand = expansion(0,0)) +
  theme_void() + # add more details later
  theme(
    plot.title = element_text(size = 9, hjust = .5),
    legend.position = "top",
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8),
    legend.box.margin = margin(0,0,-10,0),
    plot.margin = margin(0,5,0,0, "pt"),
    text = element_text(family = "Times New Roman")
  )

cities <- data.table(
  label = c("Zurich", "Geneva", "Basel", "Lugano"),
  xend = c(683240, 499950, 611260, 717875),
  yend = c(247770, 117975, 267260, 96900),
  x = c(775000, 512000, 525000, 780000),
  y = c(247770, 84000, 267260, 96900),
  curvature = c(-.2,-.4,-.3,-.1),
  nudge_x = c(1000,1000,0,1000),
  nudge_y = c(000,0,0,0),
  hjust = c(0,0,.5,0),
  vjust = c(0.5, .45, 1, .35)
)

annotater <- function(data){
  out1 <- geom_curve(
    data = data,
    aes(x = x, xend = xend, y = y, yend = yend),
    curvature = data[, curvature],
    size = 0.2, arrow = arrow(length = unit(0.01, "npc"))
  ) 
  out2 <- geom_text(
    data = data[1],
    aes(x = x, y = y, label = label, hjust = hjust, vjust = vjust),
    nudge_y = data[1, nudge_y], size = 3,
    nudge_x = data[1, nudge_x], family = "Times New Roman"
  )
  
  return(list(out1, out2))
}

mapA <- map +
  annotater(cities[1])[[1]] + annotater(cities[1])[[2]] +
  annotater(cities[2])[[1]] + annotater(cities[2])[[2]] +
  annotater(cities[3])[[1]] + annotater(cities[3])[[2]] +
  annotater(cities[4])[[1]] + annotater(cities[4])[[2]]


# Plotting dfA Change -----------------------------------------------------
myColorsB <- RColorBrewer::brewer.pal(9, "Purples")[c(5,7,9)]

mapdfA <- ggplot(data = muni_geo) +
  geom_raster(data = relief, inherit.aes = F, aes(x = x, y = y, alpha = value)) +
  scale_alpha(name = "", range = c(.6, .0), guide = F) +
  geom_sf(data = country_geo, color = "black", size = 0.05, fill = "transparent") +
  geom_sf(mapping = aes(fill = dfAChangeQuantile, color = factor(is.na(treatmentGroups))), size = .025) +
  scale_color_manual(values = c("white", "transparent"), guide = F) +
  scale_fill_manual(
    values = myColorsB,
    na.translate = F, # doesn't plot the NA or put it in the legend
    name = "Change in immigrant workers (per 100 Swiss)",
    labels = c("< 1", "1 - 10", "> 10"),
    guide = guide_legend(
      override.aes = list(color = "white"),
      direction = "horizontal",
      keyheight = unit(5, units = "pt"),
      keywidth = unit(50, units = "pt"),
      title.position = 'top',
      title.hjust = 0.5,
      label.hjust = .5,
      nrow = 1,
      byrow = F,
      label.position = "bottom")
  ) +
  geom_sf(data = lake_geo, fill = "#D6F1FF", color = "transparent") +
  scale_x_continuous(expand = expansion(0,0)) +
  labs(x = NULL, y = NULL) +
  theme_void() +
  theme(
    plot.title = element_text(size = 9, hjust = .5),
    legend.position = "top",
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8),
    legend.box.margin = margin(0,0,-10,0),
    plot.margin = margin(0,0,0,5, "pt"),
    text = element_text(family = "Times New Roman")
  ) 
mapdfA <- mapdfA +
  annotater(cities[1])[[1]] + annotater(cities[1])[[2]] +
  annotater(cities[2])[[1]] + annotater(cities[2])[[2]] +
  annotater(cities[3])[[1]] + annotater(cities[3])[[2]] +
  annotater(cities[4])[[1]] + annotater(cities[4])[[2]]

figure1 <- mapA + mapdfA + plot_layout(nrow = 1) + plot_annotation(theme = theme(plot.margin = margin()))


ggsave(filename = here("figures", "Fig1.png"), plot = figure1, width = 6.26894/1.2, height = 6.26894/1.2*0.42, units = "in", dpi = 600)
